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Figure 1. Computational Cell and h Matrix 

The position vector of atom i, rj, is in the Cartesian basis 

5=xi£i+yje2+zj§3. [1] 

The underlines denote vectors. Figure 1 shows the molecular dynamics computational cell. The 
position vector of atom i can be written in terms of the a, b> and £ vectors spanning the 

1 

computational cell edges as 

3=§3+iTife+0c |2] 

where t|, £ e [0,1] for all atoms in the computational cell. This representation is useful in 
applying penodic boundary conditions, where an atom leaving the computational cell on the 
right is considered to reenter the cell on the left. Scaling, the change of basis from the Cartesian 
orthogonal unit vectors to the a, b. £ basis is accomplished through 

-1 

3i = h ij, [3] 

and the reverse transformation through 

Q = hg . [4] 

The h matrix defines the metric tensor G through 

G = h T h [5] 

and G modifies the dot product through 

g*sj = gG^. [6] 

Also, G defines the strain tensor through 

! T-1 -1 

£ = 2 K h 0) Gh 0 ]. [7] 
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where hg is the h matrix for the reference state. 

Newton's equations of motion are, in the Cartesian basis, 



Integration of this set of coupled ordinary differential equations, the molecular dynamics method, 
provides the position of each atom as a function of time. In the scaled basis, 


$ 

§i = ^ *G G$, [9] 

and an additional matrix equation tracks changes in the size and shape of the computational 
cell through 


h = P)o - hZ] , [10] 

where W is a pseudo-mass assigned to each element of the h matrix, p is related to the internal 
stress tensor, P is applied pressure, S is related to applied stress, and s is the inverse of the 
transpose of h. 


The fluctuation of two stress or strain tensor components is just their covariance' 
5(xy) s <{x-oo)(y-<y>)> = <xy> -<xxy> , [11] 

where < > denotes a time averaged quantity. The usual formula for the mean, 


1 N 

<o = R Xxfc, 
k=1 


[ 12 ] 


risks inaccuracy through cancelation when the sum becomes large. A better formula, 

k-1 1 

<x>k = -£-<x>|<-i + £Xk [13] 


updates the average at each timestep k. 

2 

The elastic constants of a material are the tensor components C yy, the stiffnesses, or the 
components of the inverse tensor, Sp, the compliances. The relation 


Sykl - |^jT^ e li £ kl) I 14 ) 

provides values of the compliances from the fluctuations of the strain components, the reference 
state volume, and the temperature. Stress fluctuations are related to the elastic constants by a 
similar formula. 

Calculation of elastic constants requires initializing a lattice of atoms, an h matrix defining the 
computational cell boundaries, and the parameters defining the interatomic potentials, then 
numerically integrating the equations of motion for the atoms and the computational cell 
boundaries, calculating strains and stresses at each timestep, and updating at each timestep the 
average quantities required for calculation of fluctuations. 
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Applications of interest to NASA include calculation of elastic constants near lattice defects 
such as dislocations and grain boundaries, and investigation of the effect of changes in elastic 

3 

constants, caused by core ionization in XPS, on asymmetries in spectral peaks . 
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Enhancements of the DYNAMO FORTRAN code included addition of documentation, 
improved memory management and data flow, and enabling use of several interatomic 
potentials. 

Documentation included headers for each subroutine which document the parameters input 
to the routine, algorithms employed, and output. Numerous comments were added to the body 
of each routine as well. 

Improvements to memory management consisted of removing all common blocks, replacing 
them with dimension statements, thus eliminating many variables and arrays from the 
specification section of subroutines in which they were not referenced. Passing all parameters 
explicitly through argument lists clarified data flow. 

Enabling the use of all available interatomic potentials required inclusion of data structures 
for each type of potential. Setting an input variable selects the potential, obviating maintenance 
of a separate version of the code for each potential. 
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